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1 INTRODUCTION 



This is the second in a series of two papers devoted to a 
relativistic computation of torques from an external per- 
turber on a thin disc due to interactions at the Lind- 
blad resonances, i.e. locations in the disc where the orbital 
frequency Q and the radial epicyclic frequency k, satisfy 
hi — ±m(fi — fi s ), where Q s is the pattern speed of the per- 
turbation. Such resonances h ave been extensively studied in 



the nonrelativistic case (e.g . iLvnden-Bell fc Kalnajsl Il972t 



Goldreich fc Tremainelll97Sl . 1 19791 . Il98d : I Lin fc Papaloizoul 
19791 ). In the first paper ("Paper I"), we performed this com- 
putation for a general time-stationary, axisymmetric, space- 
time with an equatorial plane of symmetry and a metric per- 
turbation h a p that respects the equatorial symmetry. This 
paper ("Paper II") completes the evaluation of the Lind- 
blad torque in the case of most interest: the perturbation 
of the accretion disc surrounding a Schwarzschild or Kerr 
black hole by a small secondary also orbiting in the equa- 
torial plane. Such computations of the Lindblad resonant 
strengths may be relevant in the context of electromagnetic 
counterparts to binary black hole mergers , particularly if an 
inner disc is involved l|Chang et al.ll20ld ). (The more com- 
plicated case of perturbations outside of the equatorial plane 
- as may occur in the case of a merger where the primary 



ABSTRACT 

We present a fully relativistic computation of the torques due to Lindblad resonances 
from perturbers on circular, equatorial orbits on discs around Schwarzschild and Kerr 
black holes. The computation proceeds by establishing a relation between the Lindblad 
torques and the gravitational waveforms emitted by the perturber and a test particle 
in a slightly eccentric orbit at the radius of the Lindblad resonance. We show that our 
result reduces to the usual formula when taking the nonrelativistic limit. Discs around 
a black hole possess an m = 1 inner Lindblad resonance (ILR) with no Newtonian 
Keplerian analogue; however its strength is very weak even in the moderately rela- 
tivistic regime (r/M ~ few tens), which is in part due to the partial cancellation of 
the two leading contributions to the resonant amplitude (the gravitoelectric octupole 
and gravitomagnetic quadrupole). For equatorial orbits around Kerr black holes, we 
find that the m = 1 ILR strength is enhanced for retrograde spins and suppressed 
for prograde spins. We also find that the torque associated with the m ^ 2 ILRs is 
enhanced relative to the nonrelativistic case; the enhancement is a factor of 2 for the 
Schwarzschild hole even when the perturber is at a radius of 25M. 

Key words: accretion, accretion discs - relativistic processes - black hole physics. 



hole is rotating and the secondary is in an inclined orbit - 
is left to future work.) 

The resonant torque formula in Paper I depended on the 
geodesic properties in the unperturbed spacetime as well as 
being proportional to the square of the absolute value of the 
resonant amplitude S^ m , which was a function of the e lm ^ 
Fourier component of the metric perturbation h a p and its 
spatial derivative h a p, T . The construction of these perturba- 
tions generally depends on the solution for the Weyl tensor 
component ^4, which may be solved using a separable wave 
equation with a source given by the stress-en ergy tensor 
associated with the perturber (|Teukolskvl 1 19731 ) ; and then 
ha/s may be obtained by applyin g a second-order dif ferential 
operator to a master pot ential dChrzanowskilll975l ). which 
may be derived from ip4 (|Waldlfl978l ). Fortun ately, for our 
compu tations there is a way to circumvent the IChrzanowskil 
( 1975) procedure: Paper I showed that the particular com- 
bination of metric perturbations we require is related to 
V (m \ the power delivered to a test particle in a slightly 
eccentric orbit by the e lm< ^ component of the perturbation. 
By replacing the perturber with an equivalent gravitational 
wave source - either incoming from past null infinity in the 
case of an inner Lindblad resonance (ILR), or emerging from 
the past horizon in the case of an outer Lindblad resonance 
(OLR) - we may equate "P^" 1 ) with the power absorbed from 
the gravitational wave. However, energy is conserved on a 
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time-independent background metric, and thus "p( m > can be 
related to the interference between the equivalent gravita- 
tional wave representing the perturbation and the gravita- 
tional wave emitted by the test particle. This allows us to ex- 
press the resonant amplitude and hence the resonant torque 
in terms of the waveforms emitted by the perturber and the 
test particle (both to future null infinity and into the fu- 
ture horizon), so that standard methods to solve for ^4 are 
sufficient. 

The outline of this paper is as follows. In Section [2l we 
introduce the Kerr metric and review the associated stan- 
dard notation. Section [3] reviews the geodesies in the Kerr 
spacetime and their description with action-angle variables, 
and Section [4] describes the compuation of the perturbation 
in the Weyl scalar ^4; while both of these subjects are stan- 
dard, there are some differences in our treatment that are 
particularly suited to the problem at hand, and we make 
frequent use of intermediate results when taking the non- 
relativistic limit, so an extended discussion is warranted. 
Section [5] presents the key new theoretical result of this pa- 
per, relating the behaviour of rp4 near the horizon and at 
infinity to the resonant amplitude S^" 1 ' . We recompute the 
resonant amplitudes in the Kepler problem in Section [6l and 
then proceed to investigate the Lindblad resonances in the 
Schwarzschild problem in Section]?] Section[8]then considers 
the Lindblad resonance amplitudes associated with equato- 
rial orbits in the Kerr spacetime. We conclude in Section [9] 



The standard Newman-Penrose basis is 

I = ^ dt + -^84, + Or, 

A ( r2 + a % , a a a N \ 
n = 2E {^^ dt+ A^~ a "J' 

ia sin 9 dt + do + i esc 9 ds 

m = -= — , and 

V2 (r + ia cos 6) 

— ia sin 9 dt + do — i esc 9 8a, 



v2 (r — ia cos 9) 
Expressions involving m can be simplified if we use 



P = 



r — ia cos f 



and p ■ 



r + ia cos ( 



(3) 



(4) 



which satisfy pp — E 1 . The Weyl scalar ip4 used to describe 
the emitted gravitational waveform is 

ip4 = —C a 1 sn a fh l3 n' y rh s = — R a 0-ysn a fhP n 1 m s , (5) 

where C a p~fS is the Weyl tensor, and the equivalence to the 
component formed from the Riemann tensor R a p 1 s is due 
to the Newman-Penrose basis conditions. 

The horizons of the black hole are at radial coordinate 



1± = M ± ^M 2 -a 2 . 



(6) 



Particles very close to the horizon (r — rh+ — > + ) rotate at 
a pattern speed of the hole's angular velocity: 



fin = 



1-VT 



2Mr h+ r 2 + a 2 



2a* M 



(7) 



2 KERR METRIC AND NOTATION 
2.1 The metric and null tetrad 

We parameterize the Kerr black hole sequence with the grav- 
itational mass M and the specific angular momentum a. We 
use relativistic units where the Newtonian gravitational con- 
stant and the speed of light are equal to unity. The dimen- 
sionless angular momentum is a* = a/M. 

The Kerr metric in Boyer-Lindquist coordinates 
|Bover fc Linda uisdll967h is 



Note that for real coordinates, m = m* and p = p* , 
where * denotes the complex conjugate; however we will 
occasionally analytically continue r to complex values, in 
which case the barred quantities are not the complex con- 
jugates of the unbarred quantities: p(r,9) = p*(r*,9*) 7^ 
p*(r,9). 

Finally, we define 



K = uj(r 2 + a 2 ) — am 



and use the angular operator 



JLl 



m esc 9 + au sin 9 + n cot ( 



(8) 



(9) 



2.2 Notation in related works 



ds 



2Mr\ 2 AMar . 2 
1 ^— ] dt — — sm 9 dtd(p 



+ 



E / E 

(r 2 + a 2 ) 2 - Aa 2 sin 2 9 



sin 2 fld0 2 



-^dr 2 + Ed0 2 , 



(1) 



where A = r 2 — 2Mr + a 2 and E = r 2 + a 2 cos 2 9. The 
contravariant metric coefficients are 

tt _ (r 2 + a 2 ) 2 - Aa 2 sin 2 9 



9 



2aMr 



AE 



A — a sin 2 



AEsin 2 # ' 



g = and g = E . 



(2) 



Our notation appears to be common in the literature but 
other examples can be found. 

• We are consistent with the metric and (where applica- 
bl e) null tetrad used in the standar d general relativity text 
bv lWaldl (|l984T ). lMisner et al.l || 1973b use "p 2 " to denote our 
E, and do not fix a normalization for the principal null vec- 
tors. 

• IChandrasekharl (|l992h uses the H signature, and 

uses "p 2 " to denote our E; "93" to denote our 0; "E 2 " to 
denote our (r 2 + a 2 ) 2 — Aa 2 sin 2 9\ and "5" to de note our 
sin 2 9. For the perturbations. Ichand rasckhar] (|l992h denotes 
the frequency by —0+ , and uses the opposite sign of K. 
Additionally, our p and p are denoted by — p _1 * and — p _1 , 
respectively. However, the null tetrad and the operators ]L n 
and Kit are the same. 
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3 TIMELIKE GEODESICS IN KERR 

We utilize the Hamiltonian formulation of the equations of 
motion for a particle. As is well-known, the action for a par- 
ticle of mass jj, is S — —fj, J dr, where r is the proper time 
along the particle trajectory. For our purposes, the fastest 
route to the torque formula is not to use the covariant repre- 
sentation of the action but rather to explicitly parameterize 
the particle's trajectory using the coordinate time t, which 
is always possible outside the outer horizon. This method, 
which explicitly keeps only the 3 physical degrees of freedom, 
is best suited to a perturbation analysis. 

As in Newtonian perturbation theory analyses, it is 
most convenient to work with action-angle variables, us- 
ing the 3+1 version of the Hamiltonian that retains no 
gauge freedom assoc i ated with the particle trajectory. 
iHinderer fc Flanagan! {2008) constructed a set of action- 
angle variables in which the particle's trajectory is parame- 
terized by proper time r, and t is promoted to a dynamical 
variable (with conjugate momentum p t = —p£)- Their ac- 
tions (J r , Je, J<p) are equal to ours, since the momenta are 
the same, however the angle variables are different since 
ours advance at uniform rate with respect to coordinate 
time and theirs advance at a uniform rate with respect to 
proper time. Thus the Fourier decompositions are also dif- 
ferent. Other works that have constructed the Hamiltonian 
for geodesic motion in 4-dimensional space have projected 
the m otion into the 3 physical degrees of freedom ()Schmidtl 
2002), but appear not to have constructed the full trans- 
formation from action-angle variables to the familiar spatial 
coor dinates and momenta, which we will need to complete 
here. iFlanagan fc Hindered |2010) considered resonances in 
inspiralling black hole binaries, but parameteri ze their tra- 
jectory in terms of the "Mino time" A = Jdr/E ljMinoll2003l : 
iDrasco et al.|[2005h . This again means that they have an ad- 
ditional conjugate variable pair not present in our treatment, 
and that their angle variables advance at a constant rate as 
measured by A rather than by t. 



3.1 Hamiltonian and constants of the motion 

The trajectory of a massive particle can be followed by pa- 
rameterizing the trajectory x l (t) where x r G {r, 9,<f>}, and 
using the action S = — t. This results in the Hamiltonian 
H — —pt , where p t is determined from the p; via the mass- 
shell condition (Paper I): 



H{t,x\pi) = 



g tl p t - sj (g u Pi) 2 - g tt g i ip i p ] - p?g u 



(10) 



The timelike geodesies in the Kerr metric are char- 
acterized by three constants: the energy per unit mass 
£ = —pt/n = —Ut', the angular momentum around the sym- 
metry axis per unit mass, C — Ptp/fJ- — U4,; and the Carter 
constant, 



Q = 2E(u ■ {)(« • n) - r 2 
which may also be expressed using 



(£ - a£y 



1C = Q+(C- OS) 2 



(11) 



(12) 



Given the three constants of the motion {£, Q,C}, it 
is possible to obtain the momenta when the particle passes 
through any spatial position (r, 9, <f>). Specifically, we always 



have eastward momentum t u = C. The southward momen- 
tum given by Eq. (7.164) of IChandrasekharl l| 19921 ). 



Ug = Q — a 2 (l — £ 2 ) cos 2 I 



C 2 cot 2 1 



(13) 



and th e radial momentum by Eq. (7.160) of IChandrasekharl 
<|l992t) . 



3.2 



A 2 u 2 = [(r 2 + a 2 )£ - aCf - A(r 2 + K). 



Actions in terms of the energy, Carter 
constant, and angular momentum 



(14) 



It is useful in integrable problems to define the action-angle 
variables. We begin by considering the actions correspond- 
ing to the r, 0, and <fr loops around the invariant torus cor- 
responding to a set of constants {£, Q,C}. The 0-direction 
is the easiest: the action is 



■ho 



1 

2^ 



P4, &(f> = P4> — [iL. 



(15) 



We define the notation = J^/p, — C. 

For the 0-direction, we use Eq. (Tl3 
loop in the (9, wg)-plane. Its area, 



ue dO, 



which defines a 



(16) 



involves an elliptic function, which however is most easily 
evaluated by numerical integration. It is convenient to switch 
to the variable z — cos 9, in which case we find 

(l-z 2 )u 2 g = Q{l-z 2 )-a 2 {l-£ 2 )z 2 {l- z 2 )-C 2 z 2 . (17) 

The turning points are found at the zeroes of the right-hand 
side, which is quadratic in z 2 . These zeroes are z 2 = z±; 
inspection of the sign of the right-hand side at z 2 6 {0, 1, 00} 
shows that the zeroes have the ordering < z'L < 1 < z\. 
These zeroes can then be found by bisection. 

We may then change variables from 9 to z; noting that 
U z = (1 — z 2 )^ 1 ^ 2 u e , we find 



Q(i -z 2 y 



The action is then 
- Q 1/2 



Je 



2tt 



1-* 

1 

z 



i-£l 

z 2 



dz 
1-z 2 



(18) 



(19) 



We solve this integral with the substitution z = 
z_ sin(i7rtanh£), where the full integral is given by 4 times 
the integral J °° d£. Written in terms of £, the integrand is 
smooth, even, and decays exponentially at large f. Summa- 
tion of the integrand in £ at points (n + |)A£ thus enables 
evaluation of the integral with exponentially small error as 
A£ — > + and iVA£ — > 00 (where TV is the number of points). 

For the r-direction, Eq. (|14[) defines a loop in the (r, u r )- 
plane, and one may again find the area 



J r = — <t> u r dr. 
2n 



(20) 



A practical solution for J r is to find the turning points r_ 
and r+ by solving the quartic equation A 2 u 2 = for r, 
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Eq. (fl4[l U Then a substitution of the form 



r+ 



■ tanh /3 



(21) 



enables one to turn the integral into one over — oo < (3 < oo 
(multiplied by 2 to get the inward leg of the trajectory), 
where the integrand is analytic in the vicinity of the real /3- 
axis and declines exponentially as /3 — > ±oo; it may thus be 
evaluated by the simple method of summing the integrand 
at equally spaced abscissae (3. 

A problem one may encounter is that there is only a 
finite range of energies [£ m i n (£, Q), £ ma x(£, Q)] over which 
bound orbits can exist. The sign pattern of the extrema can 
be used to distinguish the £ < £ m - m versus £ > £ max cases. 



3.3 Geodesic properties 

For a given value of the actions (J r , Jg, J^>), one may obtain 
the constants of the motion {£, Q, £} by inverting the equa- 
tion for the actions in terms of the constants of the motion. 
The determination of C = is trivial. The determination 
of £ and Q is harder, requiring the solution of a nonlin- 
ear system of two equations; we solve these iteratively by 
first writing a function to obtain £ (J r , Q,,J$) by bisection 
solution of J r (£ , Q, C) = J r ; and then writing a function 
to adjust Q (again by a bisection search) until we find the 
desired Je. 

We will often need the 3x3 matrix of partial derivatives 



/ OE HQ oc 
/ sX dJr- r71 



M 



8Jr 



BQ 
dJ r 

<)Q 

eg 

0J o 



\ 3J& BJ4, 3J& 



d.] r 
9-U 



(22) 



The last column of M is simply (0, 0, 1) T . The first column 
is notable for being the vector of fundamental angular fre- 
quencies corresponding to the r, 8, and (j> directions on the 
torus, (Q r , Q9, ^) T - 

It is possible to obtain M by numerical differentiation, 
but it is more accurate to obtain its inverse M _1 by differen- 
tiating the actions with respect to (£, Q, C). The last column 
(the vector of partial derivatives of J$) is simply (0, 0, 1) T . 
The second column (the vector of partial derivatives of Je) 
can be obtained using the relation 



dJe 
dA 



1 

2^ 



du z 
~dA 



Az = — 
4n 



dA 



(23) 



where A £ {£, Q, £}, and we have used the fact that u z — > 
at the turning points to set to zero terms associated with 
changes in z m m,max. The explicit expressions are 



Q(l - z 2 ) - a 2 (l - £ 2 )z 2 (l - z 2 ) - C 2 z 



2 2] 1/2 



(24) 



1 We solve the equation by first finding the inflection points (via 
a quadratic equation) and then using the bisection method to find 
the extrema. Finally a further bisection gives the roots. The sign 
pattern of the extrema determines whether there are 1 or 3 roots 
outside the outer horizon; stable bound orbits require 3 roots. 



with derivatives 



d{ul) 

d£ 
d{ul) 

dQ 
d{ul) 

dC 



2a 2 £z 2 



and 



2Cz 2 



(1-z 2 



(25) 



Near the turning points or for low inclinations, u z becomes 
small, which is an issue since it is in the denominator of 
Eq. (|23|l . We thus set z = z_tanha, perform the integral 
for < a < 00, and then multiply by 4 to get the whole 
cycle; using Eq. (|18[1 this gives 



dz 
u z 



1 



Qi/2 



v 1 - 



: sech a da. 



(26) 



For large inclinations, z-fQ 1 ^ 2 may be obtained directly; 
for small inclinations (z~ < 10~ 8 ), the equatorial limit may 
be used, 



lim 



1 



S^o+ Q}l 2 ^JC 2 +a 2 {l-£i 



(27) 



A similar approach works for the derivatives of the ra- 
dial action. In this case, we need 



8Jr 

dA 



4tt 



d(u 2 ) 



OA 



dr 

U r 



(28) 



This time, the desired substitution is Eq. 1)21 p. with which 
we find 



dr 



sech 2 /3 d/3, 



(29) 



where P(r) is the polynomial on the right-hand side of 
Eq. (|14|) . If we factor the polynomial as 

P(r) = -(1 - £ 2 )(r - r'_)(r - r' + ){r - r_)(r - r+), (30) 

where r± and r'± are the four roots then we may simplify 
this to 



dr 



u r ^(l-£*)(r-r>_)( 
The derivatives are: 



: sech p d/3. 



(31) 



d{u 2 
d£ 

d(u 2 ) 
dQ 

d(u 2 ) 
dC 



2[{r 2 + a 2 )£- aC] {r 2 + a 2 ) 2a{£ - a£) 



1 



and 



A' 

-2a[(r 2 + a 2 )£ - aC] _ 2(£ - a£) 
A 2 A 



(32) 



3.4 Particle position and momentum in terms of 
the action-angle variables 

In perturbation theory it is critical to be able to obtain 
the particle's phase space location (x l ,pi) in terms of the 
action-angle variables (Ji,^ 1 ). The generic procedure to do 

2 These are all real in the case of stable orbits since P(r) is neg- 
ative at r = 0, positive at the outer horizon r = rh+, and then 
has 3 roots outside the outer horizon. 
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this is as follows. First, for a given {Ji}, we identify the con- 
stants of the motion {£, Q,jC} on the corresponding torus. 
These three actions mutually commute: { Ji, Jj}p = 0, where 
{, }p denotes the Poisson bracket. Second, we must construct 
the angle variables. For actual numerical computation, the 
method of choice is to use the direct conditions to construct 
the mapping of (Ji,^ 1 ) — > (x r ,pi), which will depend on the 
(unknown) origin of the angle coordinates ip — (0, 0, 0) on 
each torus; and we will find a valid origin b y inspection. 

W e first use the direct conditions (e.g. iGoldstein et al.l 
2002, Eq. 9.48) to write a system of differential equations 
for x r and pi as functions of the angles for fixed J: 



&if>i 



dJi 

dui 



and 



dui 



8Jn 



(33) 



These equations can be re-written in terms of derivatives of 
constants of the motion, 



dx i 



= j2 t M 

h Ae{£,Q,c} 



\A,3 



OA 
dui 



(34) 



and similarly for Ui. These equations define a solution for if), 
except that we must choose an origin if> = on each torus; 
thus all possible solutions differ by a transformation of the 
form if) 1 -> if> i + f(J). 

Our next step is to determine an appropriate choice 
of origin, i.e. the 3-dimensional submanifold of phase space 
corresponding to tp — 0. All valid choices of angle variables 
correspond to some origin (and are related to each other 
by simple phase-shifts of the angle variables on each torus), 
but in multiple d imensions not all origins correspond to valid 
angle variables Fl l Arnold! (|l978l . §50C) shows that a (locally) 
valid choice of origin is Q 1 =constant, where (Q r ,Pi) are a 
set of canonical coordinates^ We could thus choose a par- 
ticular value of (r, 9, cj>) as our origin; but this would not 
be applicable to all orbits since there is no value of r that 
all orbits cross. We prefer to choose fixed (p r ,9,<f>), which 
is also valid since Hamiltonian mechanics does not distin- 
guish between the position and momentum variable^]; we 
take p r = 0, 9 = n/2, and cj> — 0. 

It is then necessary only to apply certain inequalities 
so that each torus intersects the if} — manifold once and 
the angle coordinates are defined globally on each torus; we 
take dH/dr < and pe < 0. This corresponds to the point 
of pericentre and ascending node at zero longitude, i.e. 



(r-,~,o) and Ui = (o,~Q 1/2 ,c) 



(35) 



Starting from ip — (0,0,0), we may use Eq. (|34[) to 
evolve the particle to any chosen angle coordinates. Since 
the construction of the torus integrates over no more than 1 

3 A trivial way to see this is to note that the direct conditions 
show a transformation ip 1 tf) % + f l (J) to be canonical if and 
only if the 3x3 matrix df l /dJj is symmetric, i.e. if / is derivable 
from a potential: f % (J ) = d^/dJj fo r some <3?(.7). 

The construction in lArnoldl lll978l) technically shows that the 
generating function for the transformation (Q 1 , Pi) — > (ip 1 , Ji) 
vanishes at the chosen origin; but inspection shows that tft = 
there as well. 

5 This argument is equivalent to applying first a canonica l trans- 
forma tion Q r = p r , P r = —r, and then the construction in lArnoidl 
lll978l l. 



cycle, even a simple integrator is sufficient (we use the 4th 
order explicit Runge-Kutta method) . 

We finally need the formulas for the partial derivatives 
of £ , Q, and C with respect to {x 1 ,Pi). For £ , this is simple: 
the partial derivatives represent the Hamiltonian flow, 



d£_ _ dH_ _ .i _ 7/ 

dpi dpi u* 



(36) 



where ut is determined from the normalization g a ^u a Uf) = 
— 1 and u a is obtained by raising indices. The derivatives 
d£/dx % can be determined from the conserved quantities, 
e.g. by taking the t-derivative of Eq. (|I3|I . 

2u e u g = [2a 2 (I-£ 2 )cos6>sin# + 2£ 2 cot6>csc 2 6i]0; (37) 

using that 8 = u e /it* = ite/(E?t'), we find 

a 2 (l -£ 2 )cos0sin6l + £ 2 cot6>csc 2fl 



<■>£ 



We also know trivially that 

-£=0. 



(38) 



(39) 



Finally, taking ^ of the t-derivative of Eq. (|I4[1 gives 

A 2 u r u r + 2(r - M)Au 2 r r = 2[(r 2 + a)£ - aC]r£r - rAr 

-{r - M)(r 2 + K)r. (40) 
One then uses r = u r '/it* = Au,-/(Eit*) to obtain: 

_ d£_ _ . 2[(r 2 + a 2 )£ - aC]r£ r 
~~ dr ~ Ur ~ AEm* ~ Em* 

(r-M)(r 2 + K) 2(r-M) a 

AEu* Em* r - 1 > 

We may find the derivatives of Q by taking the differ- 
ential of Eq. (fT3|) : 



dQ = 2 [ue dug — a 2 {l — £ 2 ) cos 9 sin 9 d9 — a 2 £ cos 2 8 d£ 



+£. cot 2 8 d£. - C 2 cot 8 esc 2 8 d8] 



(42) 

Recalling that C = u<f>, and using the aforementioned rules 
to obtain the partial derivatives of £ , we may find dQ/dx l 
and dQ/dpi. 

This argument allows us to take any action-angle vari- 
ables (t/) 1 , Ji) and construct the usual coordinates (x l ,pi). 
We have not implemented an inverse function (x l ,pi) — > 
(if)' 1 , Ji) since it is not required for this work, although we do 
not expect it to present any special difficulty. 



4 GRAVITATIONAL PERTURBATIONS 

We next describe the solution of the equations for the Weyl 
tensor component ip4 given the particle trajectory. The ap- 
proach is to use the separability of the equations to write 

du 
2^ 
(43) 

where the radial function Re,nu(r) satisfies a homogeneous 
equation (in vacuum) or an inhomogeneous equation (in 
the present case, with source). The separated equation and 
the behav i our o f the radial solutions were considered by 
iTeukolskvl (|l973T >: we will thus describe in detail here only 
the aspects that are required for either our numerical tech- 
niques or for the treatment of the nonrelativistic limit. 



p'221 *Wr)^ x (*)e 1Cm *- wt) 
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4.1 Angular eigenfunctions 

We are interested here in the solutions of the latitude 
eigenfunctions Sg'^(0) that satisfy the eigenvalue equation 
l|Hughesll2OO0l . Appendix A) 

- ei, m Sl%{B) = dlSl%{6) + cot 6 deSt%(6) 

+ X 2 cos 2 SI* (6)-2s X cos6 S%± (6) 



2ms cos 9 + s^ 



Here Sf m (9) denotes the values at 
that 



(44) 



0; we understand 



(45) 



For gravitational wave problems using the gauge-invariant 
Weyl tensor component 1/14 one requires the s = —2 harmon- 
ics with x = ow. The vertical quantum number I begins at 
^ m in = max(|m|, \s\) by convention. The solution method is 
standard and is described in Appendix [Al 

4.2 The radial equation: homogeneous piece 

The radial equation can be written as, suppressing the in- 
dices Imus, 



4f A -f]- TO = -r, 

dr V dr 



(46) 



where T i s a source term to be described later and the po- 
tential is i|Teukolskvlll973l . Eq. 4.9) 

\ -K 2 ~A\(r-M)K „ 2 2 

V(r) = A '- h 8icjr +£- 2amuj + a w - 2. 

v ' A 

(47) 

We consider first the solution of the source-free homoge- 
neous equation subject to either the boundary condition of 
a purely ingoing gravitational wave at the horizon r = Th+, 
or a purely outgoing wave at r = 00. The matching condi- 
tion in between in the presence of sources will be considered 
next. 

It is standard to use the radial coordinate r* defined by 

dr* r 2 + a 2 
"dr ~ ~ ~A ' 
or explicitly (e.g. iHughesllioOol ) 

r h + r - r h+ r h _ r - r h - 



(48) 



r* = r + 



■An- 



+ 



In- 



(49) 



In the r* coordinate, the radial equation becomes 
(r 2 + a 2 ) 2 d 2 72- 2 rA-2(r-M)(r 2 +a 2 ) d72_ _ yn = Q 



A dr 2 



^.,2. -Z Inner solution 



A 



dr* 



(50) 



We consider the inner region first. In this region, as r* 
—00 and r — > r^+, we have 



r — rh+ w 2M exp 



(r* - r h+ ) 



(51) 



Then since A = (r — rh- ) (r — fh+), we have A tx e rr * , where 
r= VT^ = 2_VM^ (52) 



r l+ + a2 



In the last equality we have used the root equation for the 
horizon, r 2 1+ — 2Mrh+ + a 2 = 0. We further see that 



2(r - M) 
lim 2. 2 ~ r 

r->r h + r z + or 



and 



lim 



K 



= uj — mQw = ti7 



(53) 



(54) 



• » : h+ r 2 + a 2 

Then in the limit r* — > —00, the differential equation be- 
comes 

_ 2T^- + (vj 2 + 2irro)ft = 0; (55) 
dr; dr* 

the two solutions are then exponentials, 

72.i(r*) oc e (-^+ 2r )^ an d n 2 {r„) = e iror \ (56) 

iTeukolskvl \ 19731 ) obtained these solutions and found that 
72i corresponds to the ingoing wave and 72.2 to the outgo- 
ing wave. Thus, interior to any matter sources, the physical 
solution must be that which matches to ai72i, where ai is 
some (possibly complex) constant. Since T > and 72i in- 
creases exponentially outward relative to 72.2, no numerical 
difficulty arises in starting at some large negative value of 
r*, setting 

TZi = AV tor * and ^p. = + 2 r) A 2 e" iror * , (57) 
and integrating outward with a standard (RK4) integrator. 



4-2.2 Outer solution 

The radial equation in the outer region (r* — > 00) is not so 
well behaved. In this limit, we find 



(r 2 + a 2 )' 



;V 



-uj 2 + 4iojr* 1 + C(r* 



(58) 



and the radial equation becomes (keeping the leading-order 



terms m r* 



d 2 72. 
dr 2 



2_<m 

r* dr, 



+ UJ 



4icj 



71 = 0. 



(59) 



This may be turned into a quadratic equation for the WKB 
wave number k with the replacement d/dr* — > ik; the solu- 
tions, to lowest order in r*" 1 , are 



7 3i . , 1 

K3 = uj and K4 = — uj H . 

r* r* 



(60) 



This implies an imaginary logarithmic divergence of the 
phases, or equivalently a power-law behavior of the real parts 
of the solutions at r* ^ 00, 



72.3 — > r*e'" r * and 72.4 — > r* e "* 



(61) 



Here 72.3 corresponds to a purely outgoing wave and is the 
physical solution in problems where there is no incident 
gravitational radiation. (72.4 c orresponds to a pure l y in going 
wave.) However, as noted by IPress fc Teukolskvl (|l973l ), if 
one integrates from large to small r*, the 72.4 solution grows 
relative to 72.3, so it quickly begins to dominate. Several so- 
lutions to this problem exist in the literature, such as using 
a highly accurate int egrator such that the 72. 4 solution re- 
mains sub dominant (|Press fc Teukolskvl 1 19731 ); or evolving 
a linear combination of 72. and d72/dr* that eliminates the 
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subdominance of IZ4 as r — > 00 dPress fc 

or lacks the long-range imaginary part of the pote ntial that 

causes the divergence (|Sasaki Sz Nakamurdll98^ lbl). 

An alternative, which we use here, is to note that 
Eq. (|46|) is a regular linear ODE with analytic coefficients 
except at r £ {rh-, 00} • Therefore, if we desire IZ3 and 
d7?-3/dr at any real value of r > Th+, it is permissible to inte- 
grate the ODE on any convenient path through the complex 
plane. We note further that if 5Rr is large, then while \TZs\ 
grows more rapidly than \TZ4,\ on the real axis, if Or is al- 
lowed to be positive then \lZi\ is exponentially enhanced 
relative to |72.3|. This suggests that one may integrate not 
along the real axis itself but along a contour in the first quad- 
rant of the complex plane that begins at large r where an 
asymptotic solution is valid, and ends on the real axisQ For 
concreteness, we note that the ratio of solutions obtained 
from Eq. (|61|) should, for large r*, be 
1/2 



Ha 



oc Irje 1 " 1 "* I oc [(5Rr*) 2 + (9r*) 2 ] exp(— ojQt*). (62) 



To evaluate TLz at some real ro, we integrate along the path 

(63) 



3 . Kr 
^Jr = mm < — m — , !Kr 

for which IZ4 dominates as 5ftr — s> 00. In practice, we follow 
such a path directly to ro = r if we desire lZz(r) at r > 
Tf = max{|oj| _1 , 3Af}; for r < rf, we integrate first to rf and 
then leftward along the real axis. We have experimented 
with both a complex RK4 integrator and a Bulirsch-Stoer 
methocfl; we have used the Bulirsch-Stoer integrator here 
since it is slightly faster for similar accuracy, but we found 
both methods to be workable. 

The startin g point for the inte g ration is initialized in 
accordance with IPress fc Teukolskvl j 19731 . Eq. D15) using 
terms through order r -2 (i.e. C2); our default starting value 
ofKr is 1250max{M, M" 1 }. 



4.3 Source term 

We next need the source term Te m u>{r) in the Teukolsky 
equation. This is given bjQ as 



Ttmu(r) = / dtTi, m (r,t)e lut , 



(64) 



where 



Ti, m (r,t) = A 2 (r){yl ,5(r-ro)+a r [A 1 5(r-ro)] 

+d 2 r [A 2 6(r-r )}}. (65) 



6 Since r » r* in the large-radius regime, we may construct the 
path of integration in either plane. Here the r-plane is more con- 
venient because we have explicit analytic expressions for the ODE 
coefficients, so they can be found without writing a routine for 
the complex function r(rv) or expending the substantial compu- 
tational resources to evaluate such a function at each integration 
step. 

7 The implementation of the Bulirsch-Stoer method involved tak- 
ing steps of A5Rr = —0.03 min(SRr, u>~ 1 ). Each step was computed 
using the modified midpoint method with N = 4, 6, 8, and 16 
substeps, a nd extrapola t ed to N = 00 using a cubic polynomial 
in JV" 2 ; see lPress etahl lll992l , §§16.3 ,16.4). 

8 These equations arc provided by Min o et al and in 
slightly different form bv lHuehej l|2000l . Eq. 4.39). 



Here ro denotes the radial coordinate of the particle at 
time t, and the ^-coefficients are given as follows: for the 
5-function, 

A = [X\lL\S + 2iap sin9 E\S\C nn 

p A A z 



IK 2r\ 2a 2 sin 6 cos 9 



A 



j C !'l ft 



P_ 



K 2 „. K .2Aur-2K(r-M) 



SCrr 



(66) 



for the derivative of the <5-functior|f], 
2^/2 



Ai 



P 3A 



„ * „ 2a 2 sin 9 cos 9 _ 
ElS-\ S 



Cnm 



2p ( .K\ ^ 
and for the second derivative of the 5-function, 



(67) 



(68) 



where we have suppressed the arguments of the spheroidal 
harmonic S = Sg^ au (9,<j)). The coefficients of the stress- 
energy tensor are 



Cab = p- 



u ■ a)(u ■ b) 



(69) 



where a and b are null vectors (either n or m). The values 
of 1L[eIs and ]L\S can be obtained from Eqs. ()A8|) and 

Now for a quasiperiodic trajectory along the torus, we 
may write Ti, m {r,i) as a function of the angle variables, 
7i, m (r, t) — Te im [r,ip(t)], where each if>i(t) advances at the 
rate ipi — fij. Then we take the Fourier transform, 



Te,m(r,tp) = y)Ti,m(r\q)B lq " 4 ' , 



(70) 



where q is a lattice vector (i.e. q r , qg, and q$ are all integers). 
Using %p = -0^°' + Sit, we may integrate Eq. (|64p to get: 



(71) 



With Eq. (|7ip . we may evolve each value of q separately, 
treating 7i, m (r\q) as the source, and then sum the resulting 
perturbations. 

The Fourier components Te,m(r\q) may be evaluated as 
follows. We first see that 



Te, m {r\q) 



d 3 -0 
(2^3 



Ti,m(r, i/>)e 9 " 4 '. 



(72) 



Now if we increment t/j^ by some amount Sijj^, then it is easy 
to see that <f> is increased by 8tp^ while the other phase space 
coordinates {r, 9, u r , Ue, u<t,} remain fixed. Thus Te,, n (r,tp) 
is multiplied by exp(— imStp^). Since the complex exponen- 
tial in Eq. (|72[) is multiplied by exp(iq,j,5tp^), it follows that 



9 There is a spur ious factor of p in the second term of Eq. (4.40d) 
of lHughesI <200d) . 
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7^ m (r|q) is nonzero only if = m. In this case, the tjj^ 
integral is also trivial, so we find 



%, m (r\q) 



,q<t, 



dtp r dip 6 



7i, m (r,i/>)e 1 



q-lf, 



(73) 



(2tt) 2 

This provides a means of computing 7l,m(r\q) while doing 
only a double integral over the torus instead of a triple in- 
tegral. In practical computation, the integral is computed 
as a discretized sum over N r Ne equally spaced points on 
the ip^ = subtorus. This completes the approximation of 
Ti,m{r\q) by a finite sum over ^-functions and their deriva- 
tives. 



4.4 Solution to the inhomogeneous radial 
Teukolsky equation 

We solve the full radial Teukolsky equation via a Green's 
function method. The starting point is to recognize that 
given the boundary conditions, the solution must satisfy 



, Z dm Ki(r) r<n 

'<W-i Z out^ 3(r) r>ri 



(74) 



where Z down ' out are undetermined constants. We now sup- 
pose that the source contained a 5-function at some radius 
n, i.e. we had an inhomogeneous equation, 



VIZ - 



-5(r — r±). 



(75) 



dr \ dr 

This would imply the jump conditions that R be continuous 
at r\ and that its derivative jump by 

1 



R.'(n + e)-K'(ri-e) = 



A(n) 



These two conditions allow us to solve for Z dovjn ' ° ut : 
7 down _ TZs(ri) z out _ 72-1 (n) 



A(n)W 3 i(ri)' 



(76) 



(77) 



A(ri)W 3 i(n) 
where the Wronskian is 

W 3 i(r) =Tl 3 {r)'R! 1 {r)-ni{r)'R^{r). (78) 

The Wronskian of the two solutions to a second-order ODE 
may be obtained by elementary means: in this case, we have 
W(r) oc A(r), so we write W$i(r) = HA(r). An evaluation 
at one point is sufficient to determine K. 

We thus have the full solution in the interior region 
(r < r-) 

n tmu (r) = J2 Zt^Tl! (r)e iq ^ m 2-k5(u - q ■ SI) , (79) 
where integration of the Green's function gives 



zfZZ = I [A n 3 - A 1 n' 3 + a 2 tz^} e iq 



(2tt) 3 



(80) 



and Aq and TI3 axe evaluated at the particle position. (The 
Ai and A2 terms are obtained similarly using integration by 
parts to move the radial derivative from the argument of T 
to the argument of the Green's function.) A similar equation 
is valid in the exterior region r > r+ for the outgoing wave 
amplitude Zg^ q if we swap Hi <-»■ H3. 

Using Eq. (|43[) . it follows that in the interior region, 



1p4=p 2^ V Rl ( r ) S (,m ( 6 ) e e 



(81) 



and in the exterior region 



/ 4 \ rvOU 



1,(0) 



(82) 



At large radii, p 4 lZs —¥ r~ 1 e K * n *. Then, since the flux 
of gravitational waves at large radii is the time-average 
of |i/'4| 2 /(47ra; 2 ), we may integrate over the sphere (using 
J \S\ 2 sin^d^d^ = 2-k) to get the emitted power to 00: 



E 

i ,m,q 



Zout I 
, lm,q\ 

2u 2 



(83) 



where to — q ■ CI. The power emitted into the black hole was 
derived bv lTeukolskv fc Press! (|l974h ; the solution is 



= E < 

l,m,q 



Zdown I 2 
, lm,q I 

2a; 2 ' 



where 



8192M 5 rg + ix7(ix7 2 + T 2 )(vj 2 + 4r 2 

\c\ 2 



(84) 



(85) 



Here C is the Starobinsky- Teukolsky coefficient, whose 
squared absolute value is 

\C\ 2 = [(A + 2) 2 + 4mx-4 X 2 ](A 2 + 36mx-36x 2 ) 

+48(2A + 3)x(2x - m) + U4tv 2 {M 2 - a 2 ), (86) 

and we have used x ~ aoj an d vj — u — mQn E3 

The energy and angular momentum radiated (both to 
infinity and into the hole) are required in order to fol- 
low the evolution of c ircular or equatorial orbits un d er ra- 
diation reaction (e.g. iDetweiieil Il978l ; IShibatal [l99l . Il994l ; 
lKenneficld[l998l : iHughes! [2OO0I ) : comparison of £ and C to 
literature values can be used a test of our code. Evolution 
of generic orbits that are bo th eccentric and inclined would 
also require a relation fo r Q (|Minoll2003l ; iHughes et aill2005l ; 
iDrasco fe H ughes 20061) ■ which is not required for this pa- 
per. 

We have tested our code by checking our computed en- 
ergy and an gular momentum fluxes a gainst the results from 
Table VI of IDrasco fc Hughes! (|2006l ). for M = 1, a* = 0.9, 
semilatus rectum p = 2/(rZ 1 + r^ 1 ) = 6, and a range of 
eccentricities e [defined by r+/r- = (l + e)/(l — e)] and in- 
clinations flinc = 7r/2 — #mi n . We consider all modes with 
ma.x{£, \q^,\, \q g \, \q r \} < j max , and expect convergence as 
jmax — > 00. Comparisons are given in Table [T] 



5 THE RESONANT AMPLITUDE 

Having now solved for 1/14, it remains to compute the reso- 
nant amplitude S^ m ' from Paper I. While it would in princi- 
ple be possible to compute the metric pert urbation directly, 
by co nstructing the maste r potential ^ ( Waldl Il978l ; lOril 
2003) and then utilizing the IChrzanowskil (|l975l ) procedure, 
we will find it more useful to express l S' m - ) directly in terms 
of 1p4 . 



10 Note that iHuehei ||2000| , Eq. 4.18) contains a missing factor 
of m in the first term; [(A + 2) 2 + 4ao; m fe 
[(A + 2) 2 + 4mau) mk — &a 2 uS^ J. Also note that 
iTeukolskv fc Press! lll974l) is T/2 here. 



4a?ui 2 nk ] should read 
as defined in 
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Table 1. Comparison of our energy and angular momentum fluxes at oo and at the horizon to those o fD rasco fc Hughed l|2006h (DH) 
for M = 1, a* = 0.9, and semilatus rectum 6. The "Error" column gives the maximum fractional error of any of the four columns relative 
to DH. 



e 




J max 










i H /M 2 




Error 


0.1 


20° 


4 


-4.2574E- 


6 


+5.8126E- 


-4 


-6.7238E-5 


+8.4497E-3 


1.1E-2 






6 


-4.2576E- 


6 


+5.8700E- 


-4 


-6.7241E-5 


+8.5310E-3 


6.8E-4 






8 


-4.2576E- 


•6 


+5.8738E- 


-4 


-6.7241E-5 


+8.5362E-3 


3.6E-5 






DH 


-4.2576E- 


•6 


+5.8740E- 


-4 


-6.7241E-5 


+8.5365E-3 




0.3 


40° 


4 


-5.8169E- 


-6 


+7.0118E- 


-4 


-1.0006E-4 


+7.6189E-3 


3.7E-2 






6 


-5.8857E- 


-6 


+7.2361E- 


-4 


-1.0061E-4 


+7.8091E-3 


4.4E-3 






8 


-5.8882E- 


-6 


+7.2636E- 


-4 


-1.0063E-4 


+7.8316E-3 


5.8E-4 






DH 


-5.8882E- 


•6 


+7.2678E- 


-4 


-1.0063E-4 


+7.8350E-3 





Furthermore, since we are considering Lindblad reso- 
nances, the metric perturbations are required only in the 
interior and exterior regions, i.e. at radii r < r— or r > r+, 
where the vacuum Einstein equation is obeyed. This will 
simplify our task greatly. 

The key to the computation of the resonant amplitude 
is the result from Paper I that 



2i 



V 



(m) 



(87) 



where p( m ' is the power provided by the m Fourier mode of 
the metric perturbation to a test particle of mass fii — ¥ on 
an orbit that is slightly eccentric, oscillating between R — e 
and R + e, where e is small. 

The power can be computed without direct knowledge 
of the metric perturbations, but it breaks into two similar 
cases for the ILRs and OLRs. In both cases, we use the 
fact that knowledge of ip4, m a neighborhood around the 
test particle's radius enables determination of the metric 
perturbations (up to gauge modes and to the zero- frequency 
a £ = and 1 modes" corresponding to changes in the mass 
and spin of the hole, which provide no power) and hence 
the power is the same as that which would be provided by 
a pure gravitational wave solution with the same t[>4. 

The perturber in our case is on a circular equatorial 
orbit, hence J r = Jg — and no q r ,qg need be con- 
sidered. The mode of interest has q^, = m, uj = mQ s , and 
pattern speed Q s = fl^ (evaluated at the perturber posi- 
tion). Without loss of generality, we set the initial longitude 

^ = o. 



5.1 Inner Lindblad resonances 

In the case of an ILR, the Weyl tensor component ipi is given 
by Eq. (|81|) . This is exactly the same as the case of an in- 
coming gravitational wave with azimuthal quantum number 
in and frequency ui = mfl s with the specified amplitudes in 
each I mode. In such a situation, one may see that the radial 
mode is 

Hemu(r) = Zft° wn 7?.i(r) = ^,° wn [ci 3 7?.3(r) + cuH^r)], 



where C13 and C14 are constants evaluated in Appendix iBl 
The power in incoming gravitational waves, outgoing waves, 



and waves going down into the hole are given by 

(2c) 8 icu^rf 



Pout = 22 



2ui 2 



\C\ 2 



I „ rvdown 1 2 
|Cl3^ m 



2cj 2 



and 



E 



a\Zf. 



2uj 2 



(89) 



Now we consider our test particle. It too emits gravitational 
waves, including a set of modes at azimuthal quantum num- 
ber m and at the frequency 



uj = m£l(R) — K = mf2 s 



(90) 



These waves are emitted both down into the hole and out to 
infinity, with amplitudes Zy™ and Z° 'f m that are calculable 
by the same procedure as for the perturber, but this time 
with Fourier modes {q r ,qe,q4>) = (— 1,0, m). 

We may now obtain the power absorbed by the test 
particle using conservation of energy. There is a correction 
to the power escaping to 00 and down the black hole in 
accordance with 



SPout = 



SP d . 



KJ2 

tm 

kJ2 



ydown *7C>ut* 



rydovjn ydown* 



and 



(91) 



The power absorbed by the test particle is the negative of 
this, which can be found by expanding the real part as one- 
half the sum of a quantity and its complex conjugate: 



rn ) 



<5 Pout <5 ^down 
j- X ^ ydown 

1 \ 1 ydown* / * lyout . rydovjn\ 



ryO\lt* . ydowil* \ 
c 13^1,£m + aZj l,tm j 



(92) 



We may identify the individual contributions P^ m ' by not- 
ing that it is linear in the m Fourier mode of the metric 
perturbation; and thus it arises from the terms proportional 
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j-„ 7down ydown: 



El Therefore: 



V 



(m) 



1 X 1 ydown / yout* . ydown 

I 

9^2 E^ 



7down* 

J l,-m. 



* yout , ydown \ ( CiO\ 

c 13- Zj l,i,-m + a -^l,e,-m) ■ (Vd) 



Here Ci3_ refers to the coefficient for negative values of m 
and oj: ci3-((, m, w) = Ci3(£, — m, — a;), and similarly for a_ 
(note that the a-coefncients are real). Inspection of the ra- 
dial equation shows that Ci 3 _ = C13 and a- = 01. In the 
particular case where both the perturber and the test parti- 
cle are in the equatorial plane, there also exists a reflection 
symmetry of the emitted waveform across the equator, e.g. 
Zt™£* = (-l) m ZtZ a - Therefore the two terms in Eq. (93]) 
are equal. Thus we see that the power absorbed by the test 
particle in all of the frequency ui modes is 



-T->(m) 1 \ ^ ydown / yout* . ydown* 

r — — — "I'm y c 13^1,£m + aZj \,tm 



(94) 



This has the correct dependences: it is manifestly linear in 
/ii , which is essential since the computation of the resonant 
amplitude requires division by /ii , and also it is linear in the 
epicyclic oscillation amplitude e since the order q r Fourier 
mode of the gravitational wave scales as e' 9r '. 

5.2 Outer Lindblad resonances 

A related argument applies to the OLRs. This time, we con- 
sider a perturber on a circular orbit, again emitting at fre- 
quency ui — mils, and a test particle on a slightly eccentric 
orbit emitting at frequency 



= m£l(R) + k, 



(95) 



i.e. we are considering the (q r ,qe,q<f,) = (1,0, m) Fourier 
mode on its torus. This time, since we are considering a vac- 
uum solution outside the perturber's orbit, the perturber (or 
at least its m 7^ part) may be replaced by a gravitational 
wave coming out of the hole's past horizon. The radial mode 
amplitude is now 

n lmuj {r) = Z^n 3 (r) = Z f d r n [c3ifti(r) + c 32 TZ 2 (r)]. (96) 

The changes in power escaping to infinity and going 
down into the hole are now 



5Pou 



SP d , 



= ^E 

- ^E 



Zout yout* 



and 



Zout ydown* 
- (m Z l,£m 



(97) 



but we note that Eq. (|B14|) implies acz\ = — c* 3 . The power 
absorbed by the test particle from the m Fourier mode of 
the metric perturbation is now 



rr~y{m) I \ * yOUt / yOUt* * yd( 



(98) 



11 Since 1/14 is a complex quantity whose real and imaginary parts 
encode different components of the Weyl tensor, perturbations in 
the metric tensor, curvature, etc. are not linear in 1/14 alone but 
rather are linear in 1/J4 and Thus the m Fourier mode of 
the metric perturbation depends on both the m and — m Fourier 
modes of tp^. 



Equations (|94|) and (|98l) at first appear remarkable: they 
show that the torques at the Lindblad resonances, which 
depend on 5 (m) , can be related to the overlap between the 
gravitational waveforms emitted by the perturber and a test 
particle at the location of the resonance. But this could have 
been expected: the same time-dependent multipole moments 
that are responsible for the gravitational wave emission also 
generate resonant torques. 

We are now ready to compute the resonant amplitudes 
<;( m ) w e consider three cases. First we review the case of 
a Keplerian disc, showing how the Lindblad torques can be 
treated via the Teukolsky formalism. Then we consider a disc 
around a Schwarzschild black hole wit h a perturber, simila r 
to the physical situation envisag ed bv lChang etld] (|2010h : 
this is the first case for which the relativistic machinery de- 
veloped in Paper I and here is actually necessary, and we 
find an additional m — 1 ILR with no Newtonian Keplerian 
analogueQ Finally, we compute the resonance strengths in 
the case of an equatorial orbit around a Kerr black hole. 



6 RESONANCES IN THE NONRELATIVISTIC 
LIMIT 

The problem of Lindblad resonance torques in Newtonian 
Keplerian discs (i.e. discs in nonrelativistic motion around 
a central point mass with negligible pressure gradient) has 
been treated many times; here we treat it using the Teukol- 
sky equations. We wish to find |>S' m '| 2 for each resonance. 
This requires us first to find ^°^f'q° wn for both circular or- 
bits (the perturber) and slightly eccentric orbits (for the 
test particle). We work at radii 3> M. The solutions for 
the radial Teukolsky functions in this regime are described 
in Appendix [C] the angular functions are simply the spin- 
weighted spherical harmonics. As is well-known, the Lind- 
blad resonances can be found at values of the test particle 
radius 



n = 



m =F 1 



2/3 



ro 



(99) 



where the upper and lower signs refer to the inner and outer 
Lindblad resonances. 



6.1 Emitted waves: circular orbit 

We consider first a particle on a circular Keplerian orbit 
at radius ro 3> M, orbiting at angular velocity fi^, = 

1/0 3 /2 

M ' r . The required stress-energy coefficients phased 
to zero longitude are 

M -i/LtM 1/2 „ fiM 



C C - - ZlEM^l an d C- - 

4?o 2y/2 rl /2 



2rg 



The leading-order source term is then 

A = -± s E J \iLls{e = l) 



(100) 



(101) 



(the Ai and A2 terms have powers of r 3 ^ 2 and r 1 respec- 
tively; when they are integrated, the additional d r or <9j? 

12 The new ILR does however exist for any Newtonian potential 
with an ISCO. 
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makes these subdominant to Ao). We will find it convenient 
to define 



(102) 



so that A = -fiy tm /(2ro). 

Now for the circular orbit, a particular m-mode is ex- 
cited only at lj — mfi^ = mM 1/ ' 2 r ( ^ 3 '' 2 , and the Fourier 
mode of the torus that excites it is {q r ,qe,<l4>) = (0, 0, m). 
The downward and outward radiation amplitudes are ob- 
tained from Eq. (|80f) . with the formulae for 72a, 7^3 > and H 
from Appendix [Cl 



K 3 (r )A 



£m:l),0,m — 



fm:0,0,m — 



and 



N 2{2£ + l)k ir e +1 

fti(r )A _ . 2 - e (l-2y. f iy em (2ujy +2 r( ) 



2-{2i+ 1) 



(103) 



Note the tx r ' +1 ' and tx r§ radial behaviour; this is ex- 
pected for sourcing the order-£ multipole. 



6.2 Emitted waves: eccentric orbit 

We now consider a test particle of mass \i\ orbiting at radius 
ri, and with slight eccentricity e/ri such that the particle 
oscillates between t\ — e and n + e. We are now interested in 
the (=pl, 0, m) Fourier mode (where as in Paper I the upper 
sign represents the ILR and the lower sign the OLR), which 
has frequency ui = (m=p l)M 1//2 r^~ 3,/2 . As this is a resonance 
we will not distinguish between this value of uj and that for 
the perturber. 

The computation of Ao and negligibility of A\ t 2 proceed 
in an exactly analogous way to that for the circular orbit; 
the only differences are that (i) the true radius r differs from 
its mean value n; and (ii) we must now work at general 
longitude since we no longer have trivial angle integrals. We 
find 



. [X\yim -imd> 

The amplitude emitted to future null infinity is 



Zdown 
l;<miTl,0,r, 



R[i(r)A Ti ^T- ch/y 
N 6 2tt 



(104) 



(105) 



where the integrand may be evaluated at i\r = since the 
ip^ integral is trivial. The waveform emitted into the future 
horizon Z^.^i m ma y be obtained by replacing lZz(;r) with 
^i(r). 

The epicyclic motion in th e Kepler potential can b e 
found in any dynamics text (e.g. iMurrav fc Dermottfe oOO); 
expressed in our variables, it is, at ip^ = 0, 

r = n — ecosi/' r , and = 2— smtp r . (106) 

ri 

To first order in e, we then have 



" ^ \ n-l 

2. v 2 Tm r ri • 



(107) 



Therefore, we conclude that 



ydown 



e + i 



=F m ) and 



Zout _ -2 

l;ta;^l,0,m — 1 



2(2£ + l)k 1 r{+ 2 V 2 

( ,(^2)! Ml ej/, m (2o;) £ +M- 1 



2- (2£+l)! 



6.3 Resonant amplitudes 



(108) 



We are now ready to evaluate Eqs. (|94|) and (J98J) , each of 
which has two terms. We focus on the ILRs; the treatment 
of the OLRs is analogous. A comparison of the two terms 
shows that, using Eq. (|108l) and the relations in Appendix [Cl 



^down 
az, l;fm;=Fl,0,m 



\w\M 



£m;=pl,0,m 



2t+l 



< 1, 



(109) 



so the ^i^m-Ti.o.m term dominates in Eq. (I94|) . The actual 
evaluation using Eq. (|103[) as well gives 



V 



(m) 



E 



(£-2)!(£ + 2m)j/| m r < r 1 
(£ + 2)! 2^+1 



(110) 



The summation in Eq. (|110[) can be simplified using: 



(l-2)!yL 
+ 2)!(2£ + l) 



4tt 



= — / Pe(cos(f>) cos(mcj>) d(f> 



cos{rruj>) d<^ 



47r(^!) d^ j yrr 



°l/2 W 



2? cos i 



4(£!) d? £ 



5=0 
(111) 



Here the first equality arises by considering the spherical 
harmonic addition theorem, applying it to points on the 
equator at longitudes and <f), and taking the Fourier trans- 
form over 4>; the second from the generating function relation 
for the Legendre polynomials; and the third from the def- 
inition of the Laplace coefficient. With this, and using the 
Taylor expansion formula (and the fact that the Taylor se- 
ries of b^?2 begins with the order ? m term for m ^ 0), we 
find 



■p( m ) _ 



, . , 2mb\il 

°l/2 



(?m) 



(112) 



where here the ' on the Laplace coefficient denotes differen- 
tiation with respect to the argument. It follows that 



: (m) 



fl.UJ 



2m£l s Z<; m r% 



+ 2mb[^(v-)]. (113) 



The prefactor simplifies using fl s = M 1//2 /rg 2 and lj = 
mQ a , leaving us with 

5 (m) = — fabty'fa) + 2m6$ (?")]• (114) 

2r < 

This is equivalent to the result from Paper I using the New- 
tonian potential htt- 

For the OLRs, a similar argument holds: the diz^ilm* 
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term dominates over Z°i m in Eq. (|98|) , yielding 



The epicyclic frequency is 



(m) 



E 



max{m,2} 



(£-2)! (l + l + 2m)y? ro rg 
(€ + 2)! 2^+1 r f+ 2 ' 



(115) 

We then repeat the conversion of the summation to a Taylor 
series, this time using the identity b^^' 1 ) = ?6$(?) to 
relate the series in powers of ro /n to the Laplace coefficient 
at n/ro. This gives 



; (m) 



2r, 



1/2 



[-^b(7 2 "(w) + 2m6<7 2 '(<;+) 



-4(?+) 



~M/a'(°)*™l 



(116) 



where the last term arises for m = 1 because the summation 

(i) 



over modes begins at £ = 2, whereas the Taylor series of 



has a first-order term, bjy 2 '(0) = 1. This can be compared to 
the result for Paper I, where the last term was — ?„(5 m i. The 
two terms are exactly equal at resonance = 2 2 / 3 ; recall 
that the resonance is however the only location where S' m ' is 
needed. Indeed, if one does a Newtonian calculation of 
but working in the inertial frame (where the indirect term in 
the disturbing function is replaced by a term corresponding 
to the displacement of the primary), then one derives the 
last term in Eq. (|116|) in the form presented here. Of course, 
the two forms are equivalent on resonance as guaranteed by 
the gauge invariance arguments of Paper I. 



/2 



7 RESONANCES IN THE SCHWARZSCHILD 
PROBLEM 

We now come to our the first case where we explicitly com- 
pute angular momentum transport coefficients in a black 
hole spacetime: the Schwarzschild system. We first present 
the background coefficients and resonance locations, and 
then give the amplitudes. To simplify our expressions and 
avoid proliferation of "r/M", we will use units where the 
mass of the black hole is M = 1. 



7.1 Circular orbits: a review 

For circular orbits at radius r, the sp ecific angular momen- 
tum and energy of a circular orbit are (|Chandrasekhaj| 19921 . 
§19Mq) 



r r — 2 

£ = , and £ = — , 

Vr-3 ^r(r - 3) 

Their derivatives are 
r — 6 



£' = 



and £ = 



r — 6 



2 ( r _ 3)3/2 " 2r 3 / 2 (r- 3) 3 / 2 ■ 

The angular velocity is 

The conversion from proper to coordinate time is 
, di 



At 



r-3 



(117) 



(118) 



(119) 



(120) 



and the specific epicyclic impedance is 



Z = 



6 



r(r — 3) 



(121) 



(122) 



We see that the epicyclic frequency and impedance both 
vanish at the ISCO r = nsco — 6. 

We now suppose that a perturber is placed on a cir- 
cular equatorial orbit at radius r s > nsco = 6. Lindblad 
resonances of azimuthal quantum number m occur at 



D(r) = m(r 



-3/ 



0. 



(123) 



There is no simple closed-form solution to this equation. 
However, we can deduce its properties by noting that 



D'(r) 



5/2 



-m ± 



■ 8 



\Jr{r - 6) 



(124) 



Since (r-8)/s/r(r - 6) < 1, it follows that D'(r) < for all 
positive m and r > nsco- Thus we see that for each type of 
resonance (ILR or OLR) and for a given value of m, there 
is at most one solution to Eq. (|123[) . Furthermore, we easily 
see that D > for r w nsco and D < at r — oo, so 
there exists exactly one ILR and one OLR for each positive 
integer m. 

Here we note a key difference from the Newtonian Ke- 
plerian case: there exists an m = 1 ILR. Ordinarily, the 
innermost Lindblad resonance is the m — 2 ILR (mean mo- 
tion ratio 2:1), in which the test particle goes through two 
epicyclic periods in every synodic period. Due to pericentre 
precession, the Schwarzschild metric admits the m — 1 ILR, 
in which the orbital frequency of the perturber is equal to 
the pericentre precession frequency of the test particle. This 
is not a uniquely relativistic phenomenon, but can occur in 
any system whose attractive potential at small r exhibits a 
steeper than r _1 dependence, e.g. the potential in the equa- 
torial plane of an oblate planet. Indeed, there is a ringlet of 
Saturn at 1.29 Saturn radii, whose pericentre precession rate 
nearly matches the orbital frequency of T itan, and which ha s 
thus acquired a large forced eccentricity ijPorco et al.lll984h . 

7.2 Resonance strengths 

We may now compute S^ m ' by the method of Sec.[5]for each 
of the resonances. The first three ILRs are displayed in Fig- 
ure [T] where we plot the resonance location n as a function 
of the secondary location ro; and also the torque strength 
with the perturbing mass and disk density normalized out, 



N = 



T 



27rrig 2 E(ri) 



ivmw t 2 . 

: W(^T\ l 



j(m)|2 



(125) 



The normalized resonance strength as measured by iV has 
the advantage of converging to a constant in the Newtonian 
Keplerian limit, i.e. as ro — > oo, for the resonances that exist 
in this case (m 2 ILRs and all OLRs). Its departure from 
constant behaviour is indicative of relativistic effects. 

The resonance positions and strengths are tabulated in 
Table[2] The maximum value of £ used in the computation is 
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a balance between computation time and overflow avoidance 
versus accuracy. At very large £ and small uj, the determina- 
tion of e.g. H and C13 are susceptible to overflow errors due to 
the power-law behaviour with large indices (r 1_ and r 2+e ) 
of the radial solutions to the Teukolsky equation between 
r ~ 2 and r ~ cj _1 [3 Fortunately, for the results in this 
paper we do not need to work in a regime where overflow 
occurs. For most cases, have used ^ max = 20 for the compu- 
ations at 20 < r < 250 and l max = 40 at 8 < r < 200 
For the m = 1, 2, and 3 ILRs presented, we have estimated 
the truncation error in i by extrapolatin eF^l the sequence of 
contributions from successive £; such errors are found to be 
^ 0.1% (ro = 1 and m = 2) and sC 1% (m = 3). 



7.2.1 Them ^2 ILRs 

The m ^ 2 ILRs exist in the Newtonian Keplerian limit as 
(m — 1) : m mean motion resonances, and are located at a 
fixed ratio of semimajor axes, or in this case, orbital radii: 



lim 



a. 

ro 



2/3 



0.63 
0.76 



ro = 2 
m — 3. 



(126) 



These formulae would correspond in the left panel of Fig. Q] 
to straight lines with unit slope (since this is a log- log plot). 
In fact they are relatively good approximations even at mod- 
est values of ro: for the m = 2 ILR, for example, n/ro 
increases from 0.63 (ro = oo) to 0.67 (ro = 50) to 0.72 
(r u = 20). As the secondary approaches the ISCO, however, 
the resonance locations must remain between the secondary 
and the ISCO, and hence 



lim ^ = 1. 

'•o-^nsco ro 



(127) 



This behaviour can be seen in the left panel of Fig. [T] where 
all of the resonance location curves converge to the point 
(ro,ri) = (6,6). Of course, for any finite mass ratio, the 
assumptions used throughout this paper of weak perturba- 
tions and a thin disc would break down before this point is 
reached. 

The resonant strength (as measured by N) approaches 
a constant in this limit, 



lim N 

Tq—^oo 



11 < - 



-,(m),, - 

7 ' 1/2 



+ 2mb\ 



(m) (c~) 



(128) 



This evaluates to —2.36 for m — 2 and —7.50 for m = 3; 
the convergence to these constant values can be seen from 
the right panel of Fig. [T] As one moves inward toward the 
ISCO, the strength \N\ increases. The qualitative effect is 
unsurprising since the resonance locations become closer to 



In principle such errors could be removed by working with 
ln7?,(r) instead of TZ(r), but we have not done this as it would 
have resulted in much more complex code (including branching 
to avoid numerical instabilities when 1Z passes near zero). An 
alternative would have been to define a new floating-type data 
type with more bits in the exponent. 

14 The exceptions are that for 20 < ro ^ 40 we use ^ max = 40 
for the m = 3 ILRs; and for ro > 20 we use ^ max = 30 for the 
m = 2 OLRs. 

15 Since there is a strong odd-even pattern to the contributions 
from successive multipoles, we used the last two even £s to gener- 
ate a geometric sequence of even is and did a similar independent 
procedure for the odd Is. 



the secondary. It is however noteworthy that the m ^ 2 ILR 
strengths are enhanced substantially relative to the Newto- 
nian Keplerian limit even at large distances from the black 
hole: the deviation is already 10 per cent at ro = 160, and 
reaches a factor of 2 at ro = 25. 



7.2.2 Them=l ILR 

For the m = 1 resonance, the strength is however much less, 
especially in the nearly Newtonian regime. This is in part 
due to the location of the resonance, with n < ro, and also 
due to the fact that the Newtonian quadrupole tidal field 
does not contribute to S' m ': reflection symmetry across the 
equatorial plane allows only m £ { — 2, 0, 2} contributions to 
the tidal field, and so the lowest-order contribution to the 
resonance strength comes from the (gravitoelectric) octupole 
{1 = 3). 

While the m = 1 ILR does not exist in the Newto- 
nian Keplerian problem, its location and strength may be 
estimated in the large-ro limit. The m — 1 ILR location is 
determined by the condition that the pericentre precession 
rate, 



ygT-Vn -6 3 
Q-K = l ^ w — , (129) 



correspond to the secondary orbital angular velocity, r 
This implies, for large ro, 



3/2 



_ Q 2/5 3/5 _ -| rr 3/5 

n « 3 ' r ' w 1.55r . 



(130) 



One can see this behaviour in the left panel of Fig. [T] 
Eq. (|130[) predicts that the m = 1 ILR location curve should 



be a straight line with slope |, which is indeed correct at 
large ro. The deviation from this expression is only 8 per 
cent at ro = 20, which is remarkable. 

The strength of the resonance in the large-ro limit can 
be estimated from Eq. (|110p : the leading-order term is I — 3, 
which gives 



4.36r " 19/10 and N 



-19.2r, 



-13/5 



(131) 



This result is valid at very large ro. However, at even mod- 
est ro it substantially overestimates the strength of the 
m = 1 resonance: the true N is smaller by a factor of 
0.75 at ro = 250 and 0.66 at ro = 50. The principal rea- 
son is that there is another contribution to V (1 ^ from the 
gravitomagnetic quadrupole mode (£ = 2, negative parity), 
which does not exist in the Newtonian theory but has the 
correct symmetry properties for two equatorial orbits to in- 
teract via an m — 1 mode. Roughly speaking, the grav- 
itomagnetic interaction should give a contribution to 
that is suppressed by the product of the orbital velocities 



VqVi 



-1/2-1/2 



-4/5 



but (due to the angular mo- 
mentum barrier for I — 2 versus 3) enhanced relative to 
the gravitoelectric octupole by a factor of (ri/ro) _1 ~ r 2 J° . 
Thus overall, the gravitomagnetic quadrupole interaction is 
only weaker than the gravitoelectric octupole by a factor of 
~ Tq 2 ^ . It turns out that the two contributions to S 1 ^ have 
opposite sign, resulting in a suppression of the m = 1 ILR 
strength. The correction is not small: 



iV ' (magnetic quadrupole) 
<SM (electric octupole) 



-0.19, r = 250 



-0.38, 



50, 



(132) 
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Figure 1. The locations and strengths of the inner Lindblad resonances in the Schwarzschild spacctime, shown on logarithmic axes. The 
m ^ 2 resonances have Newtonian Keplerian analogues; their strength in these units approaches a constant in the Newtonian Keplerian 
regime, but grows rapidly as one approaches the ISCO. The m = 1 resonance exists only due to the pericentre precession, and is found 
at much smaller radius; it is also much weaker, although its strength grows as we move inward. 



Table 2. The resonant strengths in the Schwarzschild problem for m = 1 and m = 2 Lindblad resonances; here ro is the orbital radius 
of the perturber and r\ is the resonance location. Truncation errors due to choice of £ max are estimated to be ^ 0.1% for the cases given 
in the table. 



m = lILR(0:l) m = 2 ILR (1:2) m = 2 OLR (3:2) m = 1 OLR (2:1) 

ro n N ri N ri N ri N 



8.00 


6.48 


-7.11E-1 


6.97 


-4.69E+1 


9.55 


+2 


28E+1 


11.33 


+7.05E-2 


9.00 


6.80 


-2.32E-1 


7.57 


-2.55E+1 


10.91 


+ 1 


79E+1 


12.99 


+8.76E-2 


10.00 


7.13 


-1.10E-1 


8.19 


-1.74E+1 


12.26 


+ 1 


51E^ 


-1 


14.62 


+1.01E-1 


12.00 


7.78 


-4.08E-2 


9.44 


-1.09E+1 


14.92 


+ 1 


30E^ 


-1 


17.85 


+1.19E-1 


14.00 


8.40 


-2.10E-2 


10.70 


-8.27E+0 


17.57 


+ 1 


18EH 


-1 


21.07 


+1.31E-1 


16.00 


9.01 


-1.28E-2 


11.97 


-6.85E+0 


20.21 


+ 1 


11E- 


-1 


24.26 


+1.40E-1 


18.00 


9.58 


-8.57E-3 


13.23 


-5.96E+0 


22.84 


+ 1 


07EH 


-1 


27.46 


+1.46E-1 


20.00 


10.14 


-6.12E-3 


14.50 


-5.37E+0 


25.47 


+ 1 


04E^ 


-1 


30.65 


+1.51E-1 


30.00 


12.67 


-1.89E-3 


20.82 


-4.00E+0 


38.60 


+9 


60E^ 


-0 


46.56 


+1.66E-1 


40.00 


14.89 


-8.69E-4 


27.13 


-3.49E+0 


51.72 


+9 


30E^ 


-0 


62.45 


+1.73E-1 


50.00 


16.91 


-4.86E-4 


33.44 


-3.22E+0 


64.83 


+9 


13E^ 


-0 


78.33 


+1.77E-1 


75.00 


21.36 


-1.72E-4 


49.21 


-2.89E+0 


97.59 


+8 


94E^ 


-0 


118.03 


+1.82E-1 


100.00 


25.25 


-8.33E-5 


64.97 


-2.75E+0 


130.36 


+8 


85E-I 


-0 


157.72 


+1.85E-1 


150.00 


32.00 


-3.01E-5 


96.47 


-2.61E+0 


195.88 


+8 


77EH 


-0 


237.10 


+1.87E-1 


200.00 


37.91 


-1.46E-5 


127.98 


-2.54E+0 


261.40 


+8 


73EH 


-0 


316.47 


+1.88E-1 


250.00 


43.25 


-8.35E-6 


159.48 


-2.50E+0 


326.92 


+8 


70E-I 


-0 


395.84 


+1.89E-1 



and then the resonant torque depends on the square of 5 (m) 
so these corrections are effectively doubled. 

The reason for the opposite sign of the gravitomagnetic 
quadrupole contribution can be understood from linearized 
gravity arguments. To lowest order, a moving particle in the 
vicinity of a mo ving perturb er experiences a gravitomagnetic 
"acceleration" (jWaldll 19841 . §4.4a): 

a = -4v x B, B(r) = qMv x , V - V ° , (133) 

r — ro\ 

i.e. B is the field generated from the momentum in the same 
way that a magnetic field is generated by electric current. 



Here ro is the position of the perturber and vo is its ve- 
locity. The test particle experiences an inward gravitomag- 
netic acceleration that is strongest at inferior conjunction 
(i.e. when the longitudes of the test particle and perturber 
are equal). This is the opposite of the Newtonian gravito- 
electric octupole field, which produces an outward force at 
inferior conjunction. 

7.2.3 The OLRs 

The outer Lindblad resonances, being external to the per- 
turber, are more similar to their Newtonian counterparts 
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than the inner Lindblad resonances. The limiting strengths 
as ro — > oo for the m = 1 (2:1) and m = 2 (3:2) OLRs 
are N = 0.19 and N = 8.62 respectively; their behaviour at 
smaller radii is shown in Table [2] 

For the strong m = 2 OLR, the resonant strength in- 
creases as we move inward because the Lindblad resonances 
are closer to the perturber than they are in the Newtonian 
Keplerian case. However, the weaker m = 1 OLR (2:1) suf- 
fers from the same partial cancellation of gravitoelectric oc- 
tupole and gravitomagnetic quadrupole contributions as the 
m = 1 ILR. Therefore at small radii it actually becomes 
weaker. 



8 RESONANCES IN THE KERR PROBLEM 

We may now move on to the resonances associated with the 
circular, equatorial orbits in the Kerr spacetime. Again, we 
use units where the mass of the primary hole is M = 1, 
and hence a = a*. We consider orbits with <f> > 0; thus 
a > (prograde spin) refers to the case where the disc orbit 
and black hole spin are in the same direction, and a < 
(retrograde spin) refers to the opposite case. The machinery 
we have developed in the previous sections is completely 
general and may be used to compute resonance strengths in 
Kerr with no new difficulties. 

The problem is very similar to that of the Schwarzschild 
spacetime: there exists an ISCO at which k — > 0, and hence 
once again there exists an m = 1 ILR. This time the basic 
frequencies are 

^ 1 , ^ /, 6 8a 3a 2 ,„„.-> 

= -rmT-a ^ K = Q V 1 ~ r + ~ — (134) 

|Okazaki et alj|l987l . Appendix). The sign of the a term in 
k/Q implies that pericentre precession is enhanced for a < 0; 
the same effect is responsible for the larger value of risco 
for retrograde spin. 

In Fig. [2] we explore the location and strength of m — 1 
ILR as a function of the secondary (perturber) location ro 
and the spin of the primary a. We would intuitively expect 
that retrograde spin (a < 0) would both move the resonance 
location n outward and increase its strength. This expecta- 
tion is confirmed numerically. Moreover, the effect is quite 
strong: even at ro = 250, a spin of |o| = 0.9 leads to a factor 
of 1.17 difference in the m = 1 ILR location depending on 
the direction of the spin (n = 39.8 for prograde, 46.4 for ret- 
rograde) and a factor of 2.4 in the strength \N\ (5.2 x 10~ 6 
for prograde, 1.3 x 10 -5 for retrograde). The difference in 
resonant strength between prograde and retrograde configu- 
rations becomes greater as ro moves inward, and at ro = 20 
and \a\ = 0.9 is more than an order of magnitude. 

At very small radii, we once again have the behaviour 
ri — > ro and \N\ — > oo as ro — >• nsco- This behaviour is 
present but not obvious in Fig. [2] because risco depends on 
a (it is larger for the retrograde configuration). 

The variation of the Lindblad resonance locations and 
strengths at fixed ro but varying a is displayed in Fig. [3] for 
ro = 50 and Fig. [4] for ro = 20. For the retrograde spins all 
of the resonances move closer to the perturber, and corre- 
spondingly they are strengthened. However, we can see that 
the effect is strongest for the m = 1 ILR, which is unsurpris- 



ing since it is closest to the hole and therefore most affected 
by spin. 



9 DISCUSSION 

The Newtonian formulae for the torque applied to a disc at 
the Lindblad resonances associated with a perturber on a 
circular equatorial orbit have been extended into the rela- 
tivistic regime. The calculation has revealed both new phys- 
ical effects, and has provided a mathematical connection be- 
tween seemingly disparate phenomena: resonant torques and 
gravitational radiation. 

At the physical level, we have learned that relativistic 
effects introduce an additional m = 1 inner Lindblad reso- 
nance at which the pericentre precession rate of the test par- 
ticle matches the pattern speed of the perturbation. This has 
no Newtonian Keplerian analogue, but in quasi-Newtonian 
language one can think of it as being due to the steepening 
of the potential. Indeed, any Newtonian potential with an 
ISCO will have this resonance. We found, however, that the 
quasi-Newtonian calculation of the resonant strength, which 
is due to the tidal octupole, is suppressed by tens of percents 
due to gravitomagnetic corrections even at ro/M > 100. In 
this sense the m = 1 ILR is a relativistic beast. 

At the mathematical level, our method of computation 
has revealed a connection between, on the one hand, angu- 
lar momentum transfer via the Lindblad resonances; and on 
the other hand, the product of the gravitational wave signals 
emitted to infinity and into the hole by the perturber and 
the test particle (assuming the latter to be in an orbit of in- 
finitesimal eccentricity). This connection arose from general 
principles: (i) the conservation of energy and angular mo- 
mentum when the contribution to both from gravitational 
waves is included; (ii) the fact that, aside from the I = and 
1 modes that do not contribute to resonant transfer, the en- 
tire perturbed spacetime structure in the vacuum regions is 
determined by the radiation degrees of freedom, described 
for Type D spacetimes by ip4,; and (iii) the ability to describe 
epicyclic motion of the test particle via Hamiltonian dynam- 
ics. This was not expected when we began the calculation, 
and we are still lacking an intuitive explanation. 

The relativistic corrections to the Lindblad resonance 
formulae - particularly the existence of the new m — 1 
ILR and the strengthening of the m 2 ILRs - may be 
important in binary black hole merger scenarios that in- 
vo lve an inner disc. T his is especially true for the proposal 
of IChang et al.l (|2010l ) , in which a bright electromagnetic 
counterpart is produced by resonant heating of this inner 
disc. A more full treatment of disc evolution including the 
new resonance as well as other Newtonian aspects of disc 
physics is beyond the scope of this paper; however, simple 
considerations sugg est that this would be a fruitful exercise. 
IChang et al.l 1120101 ) computed the inner disc evolution for a 
primary hole of mass M — 10 7 Mq and mass ratio q = 0.1, 
used Newtonian formulae for the torque, and treated the res- 
onant torques as continuously distributed in radius (which 
may be appropriate for sufficiently small |ro — ri|). They 
find that the inner disc is trunca ted at ri < 0.63rp until 
r ~ 20M (see Figs. 3 and 4 of IChang et~afl |2010| ); it is 
thus plausible that in a full treatment including the discrete 
nature of the Lindblad resonances, the strong m — 2 ILR 
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Figure 2. The locations (left) and strengths (right) of the m = 1 ILR for equatorial orbits in the Kerr spacetime. The points show 
computations using our perturbation theory code, with the symbols indicating the choice of primary spin a. For progradc orbits (a > 0) 
the resonance moves inward and become weaker, whereas for retrograde orbits (a < 0) the resonance moves outward and becomes 
stronger. 
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Figure 3. The locations and strengths of the Lindblad resonances as a function of black hole spin for a perturber in a circular orbit at 
r = 50M. 
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Figure 4. The locations and strengths of the Lindblad resonances as a function of black hole spin for a perturber in a circular orbit at 
r = 20M. 



would truncate the disc. If this is the case, then even the 
weak m = 1 ILR could be a significant contributor to reso- 
nant heating: while it is 3 orders of magnitude weaker than 
the m = 2 ILR at ro = 20, if material in the m = 2 ILR has 
been mostly cleared it is no longer obvious which resonance 
dominates the torque. This is especially true for retrograde 
configurations, where the m = 1 ILR is enhanced. While the 
distribution of values of a is presently quite uncertain, in the 
context of electromagnetic counterparts to a low-frequency 
gravitational wave detector such as the Laser Interferometer 
Space Antenna the value of a for e ach event will in man y 
cases be known to high precision fe.g. lLang fc Huehesl200r3 ). 
Due to the weakness of the m = 1 ILR, it may also be impor- 
tant to account for other weak resonances, e.g. inclination 
resonances in the case of a spinning primary; we have not 
computed the strengths of inclination resonances in this pa- 
per, but note that the techniques described here should be 
applicable to that problem. 
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APPENDIX A: SPHEROIDAL HARMONICS 



This appendix considers the solution to the angular eigen- 
mode equation, Eq. ([44]), for S^(9). 

The most convenient way to solve Eq. (|44|) is to 
write the eigenfunctions as linea r combinations of the spin- 
weighted spherical h armonics |Newman fc Penrose! 1 19661 : 
iGoldberg etal1ll967h. as has been done in previous works 
(e.g. iPress fc Teukolskvll 19731 ; lHughes|[200ol ') 



E b jem (x)Yl m (e), 



(Al) 



j — max( m| , | s | ) 



where the coefficients bji m (x) satisfy the eigenvalue equation 
dPress fc TeukolskvllT973l , §IIIa), 

cb = e! >m b, (A2) 

where b is a vector of length j max — jmin + 1 where j m i n = £ m i n 
and jmax is the highest angular momentum harmonic used 
in the finite basis set. The matrix C is real and symmet- 
ric, and is band-diagonal in the sense that C, 



if 



li- 



■f > 



2 |Press fc Teukolskvll 19731 ). In numerical compu- 



tation, we truncate at j ma x, obtain the eigenvalues £| m and 
eigenvectors b by Jacobi iteration, and compute the residual 



(A3) 



jmax is increased until e falls below some error threshold 



(usually 10 ) for all desired £. The eigenvectors are normal- 
ized using ^2 j b 2 ji m — 2-7T, which is equivalent to the usual 
normalization, 



\si%{ef de = i. 



(A4) 



The spin-weighted spherical harmonics are computed 
directly from the rotation matrices, 



Y e ° m (9) = (-ir[D(9)l 



-l) m [exp(i#L 2 



(A5) 

where L2 is the angular momentum operator around the 2- 
axis in the spin-^ representation of SO(3)0 The complex 
exponential is computed by a quadratic expansion for small 
9 (9 < 10~ 8 ), and for larger values by repeated squaring 
of the rotation matrix D(9) (each squaring doubles 9). For 
this process, we actually store D(9) — 1 where 1 is the (2£ + 
1) x (21 + 1) identity matrix; this is numerically preferable 
for small 9 to avoid exponential amplification of rounding 
errors in the squaring process. The squaring process is then 



D(20) 



2[D(0) 



[D(9) - 1] 



(A6) 



This method is slow but is stable, simple to code, and does 
not suffer from underflow occurrences (common in many 
publicly available spherical harmonics routines even at mod- 
est £). It also returns estimates of the ^-derivatives with no 
extra effort since 



dD(fl) 
(18 



iL 2 D(0). 



(A7) 



For evaluation of the source terms, we require JL\S and 
Ej\Ej\S. Given S and deS, it is easy to compute 

L\S = d e S + (-m esc 9 + x sin 9 + 2 cot 9)S. (A8) 

We further see that 

1j\1j\S = 9|S + (-2mcsc6» + 2xsin6» + 3cot6»)9eS' 



i/2 2 , 

+ (m esc ( 



+ x sin 2 9 — 2 — 2mx 



-2m esc 9 cot 9 + 4x cos 9)8. 



(A9) 



We may now use the angular Teukolsky equation for 5", 
which is a second-order ODE that expresses dgS in terms of 
5", dgS, and the eigenvalue £ . Substituting out <9f S, we find 



2(-m esc 9 + x sin 9 + cot 9)d e S 
+ [— x 2 cos 29 — 2rax + 2m 2 esc 2 t 



—6m esc 9 cot 9 
which is the form we use. 



2 + 4csc tf-flS, 



(A10) 



APPENDIX B: SCATTERING MATRIX 

Here we concern ourselves with the scattering matrix relat- 
ing the four solutions of the radial Teukolsky equation, 



tti(r) 
K 2 (r) 



C13 C14 
C23 C24 



H 3 (r) 
ft 4 (r) 



(Bl) 



where the c a b are 4 complex coefficients that we wish to 
compute. (We may also want the inverse matrix.) Our goal 

16 With the standard (Condon-Shortley) phases, iL.2 is real and 
antisymmetric. 
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here is the numerical computation of the c a b coefficients an- 
alytically from K and the parameters of the problem. 

The Wronskian of any two solutions is W a b = TZaTi'b — 
IZbR-'a and is proportional to A. In particular, the asymp- 
totic solutions give at the horizon gives 



where 



W12 = 2{iw - r) 



P = 2iMr h+ uj - 



dr 



A = 2/3A, 



2VM2 



(B2) 



(B3) 



The solutions at large radius give W34 = — 2io;A. We have 
also set Wzi = NA. 

The above Wronskian elements constrain the c a b. First, 
Eq. (|B1[) sets W12 equal to W34 times the determinant of 
the matrix of c a b, so: 



C13C24 



.P 



C14C23 = 1— ■ 



(B4) 



Second, the definition of N implies that HA = — C14VK34, so 



Cl4 



2w 



(B5) 



Further relations can be found from considering the con- 
servation of energy. For a general case with 

72(r) = 6iW x (r) + & 2 7l 2 (r) = b 3 TZ s (r) + b 4 Tl 4 (r), (B6) 

the conservation of energy (jTeukolskv fc Press! 1 19741 ) then 
provides the relation 



N 2 +a|6i| 2 



(2w) 



I64I 2 + Q2I62I 2 . 



(B7) 



Here the 62 term denotes power emerging from the past hori- 
zon, whose value is not required here. This relation may be 
evaluated for the case of 72. = 72. 1 + a 722; equating terms on 
both sides proportional to 1 and a (or a*) gives respectively 



C13 + a = , „,„ |«| 



|C| 2 



and 



C23C13 



ICI 2 



C 2 4Cl4- 



(B8) 



(B9) 



Equation (|B9[) enables us to solve for C23 in terms of the 
other coefficients; substituting into the determinant relation, 
Eq. (|B4[I , eliminates C23 and generates a linear equation for 
C24 in terms of C13 and C14: 



(2o;) 8 ci4 



Cl4 + C13 



C24 



■J 



(BIO) 



Using Eq. (|B8|I and substituting for C14 (from Eq. IB5|I sim- 
plifies this to 



C24 



iPch 



and hence 



C23 



i/3 (2w) 



■c 14 . 



(Bll) 



(B12) 



auj \C\ 2 

The programme to compute the c a b is thus as follows 



• First obtain C14 from Eq. (|B5|I and the solution for N 
from Sec. l4~4l 



• Next obtain C13 by integrating the 72i solution along the 
real axis to large r, where the 723 solution becomes domi- 
nant. By dividing by the asymptotic form for 723 (again 
keeping the first two coefficients in the expansion), obtain 
the coefficient of 723 in 72i, i.e.ci3. 

• Evaluate p and then use Eqs. (|B11|) and (|B12p to obtain 
C24 and C23- 

The inverse transformation coefficients C31, C32, C41, and 
C42 can be obtained in accordance with 



I C31 


C32 




\ c 4 i 


C42 , 





C24 — C14 
^ C23 C13 



(B13) 



we note that the substitution of the formula for the deter- 
minant in the denominator is required if this relation is used 
for numerical computation because of the very large correla- 
tion coefficient of the matrix, i.e. for some practical cases we 
have C13C24 ~ C14C23. However, for formulas involving C31 it 
is more convenient to combine this with Eq. (|B11|) to obtain 



C31 = 



£13 
a 



(B14) 



APPENDIX C: 
FREQUENCY 



RADIAL MODES AT LOW 



This appendix describes the radial modes in the nonrela- 
tivistic regime, i.e. where uj M _1 and M < r < uj^ 1 . 
This is the regime relevant for Newtonian Keplerian discs 
(Section [6}. The angular modes simply reduce to spin- 
weighted spherical harmonics with separation constant £ = 
t{l+l). 

There are infinite (logarithmically divergent in r or r — 
j"h+) phase errors in our approximations here; this does not 
concern us since the absolute phases of 72i at the horizon or 
723 at infinity cancel out of the computation. 

The solution of the radial modes in te rms of 1J1 
functi ons is described in grea ter generality by [Mano et all 
(|1996T) : see also the review bv lSasaki fc Tagoshil (|2003l . U). 
We sketch here a simplified derivation for the specialized 
case of small to, which does not require a "renormalized an- 
gular momentum parameter" and has much shorter expres- 
sions. 



CI The 72i solution 

The 72i solution (no radiation emerging from the past hori- 
zon) in this regime can be constructed by taking ui — > 0. 
With this simplification, the radial Teuk olsky equation ca n 
be reduced to a hypergeometric equation (|Mano et al.lll99rih . 
The solution is 

72 oc {-~xf- w,2 {l-x)- iT/2 2 Fi(£ + l-iT,-£-iT;3-vr;x), 

(CI) 

where 



>"h+ 



rh+ - r 

r h+ - r h - ~ 2VM 2 - a 2 



(C2) 



and r = —am/\/M 2 — a 2 . Outside the horizon we have x < 
0, and we take the branch arg(l — x) = arg(— x) = of the 
fractional powers. 
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The normalization of 72i can be obtained by taking the 



limit as r 



'h+ 



(— x — > + ). This gives 



72i ->■ A 2 e ir ' 



16(Af 2 - a 2 ) 2 e im¥, °(-x) 2 - iT/2 , 



(C3) 

where ipo = |a* + QhM ln(l — a 2 .) and we have substituted 
for Oh in order to simplify the exponent of —x. We thus see 
that 



lim72i(r) = 16(M 2 - a 2 ) 2 e lmw (-x) 2 ~' T/2 (l - a;)" 1 " 



x 2 Fi(e+l-ir,-£- 



\2-1t/2/ 

ir; 3 — ir; a;) 



/2 



(C4) 



The serie s can be made finite using the linear transformation 
formula l|Abramowitz fc Stegunll 19721 , Eq. 15.3.3): 



lim72.i(r) = 16(M 2 



2-,2 Jm V0 , 



(-x) 2 - iT/2 (l-x) 2+1 



2 + ir/2 



x 2 F 1 (2-e,3 + £;3-iT-x). 



(C5) 



To reach the Keplerian regime, we must follow this to the 
regime where —x ^> 1. Taking the highest-order (r ~ 2 ) term 
in the series, we find 



-16(M 2 



a axa eiTO », ggj r(3-ir) 



(£ + 2)! r^+l-ir) 



£+2 



.2VM 2 - a 2 , 
For our purposes, this may be written as 

Hi{r) -> fci/+ 2 , 



(C6) 



(C7) 



where using the recursion relation for the F functior0, 

r e i - 1 / 2 

1/f 2 „2^(2-l)/2 (21)! 



|fci| = [2(M 2 - a 2 )] 1 



(£ + 2)! 



(C8) 

We will not require the phase of ki; indeed, the phase is 
meaningless at the level of approximation here because in 
taking ui — > we introduce a phase error of ~ M r *!> which 
diverges as one approaches the horizon. 



C2 The 72 3 solution 

We are now interested in the solutions that asymptote to a 
purely outgoing wave at r — > oo. In this cas e, we keep cj but 
approximate M, a — > 0. iMano et al.l (1 19961 ) also provides a 
solution in this case in terms of a confluent hypergeometric 
function. They find 72 oc rf, where z — tor and f(z) satisfies 
the relation 

z 2 ^- + [z 2 -4iz- £(£+l)]f = 0. (C9) 
dz z 

As is well-known, this equation reduces to a i_Fi-type series 
upon the substitution f(z) — e ±lz g(z). Four solutions may 
be obtained this way, depending on the chosen leading power 
of z: 



Of these, 72b and 72d have the advantage of having truncat- 
ing (polynomial) if\ series; due to the nature of their os- 
cillating parts they are manifestly linearly independent and 
provide a complete basis. The highest power in r shows that 
72b (z) yields the 72 4 solution and 72d(z) yields the IZ3 solu- 
tion. The normalization is easily obtained from the highest 
term: 



72 3 (r) 



(2£)\r 1 - 



(£ - 2)\ (2ui) e + 2 



1-F1 (-2 • 



-21: -2iwr) 



(Cll) 

This is only valid in the limiting case where M — > 0; finite 
mass introduces a logarithmically divergent phase error due 
to the long-range nature of the background metric pertur- 
bation (the asymptotic expansion of dr*/dr — 1 begins with 
the order r _1 term). 

For the Newtonian Keplerian problem, we require the 
near-field solution r <ti u~ , where 



72-3 (r) = ksr 1 



with 



ka = i 



(2£)\ 



(«-2)!(2w) 



1+2 ' 



(C12) 



(C13) 



A similar result allows us to normalize IZ4,: in the near- 
field zone, 72.4(r) = k^r 1 ' 1 with 



ki — i 



: (21)! 

(e + 2) 



•(2c) 2 



(C14) 



C3 Wronskians and scattering coefficients 



The Wronksian of the 72i and 72-3 solutions is easily eval- 
uated in the Keplerian range of radii. It leads to N = 
{21 + l)fcife 3 , hence 

.2 (2£+l)\ 



:ki. 



(C15) 



(£ - 2)! (2ujY+ 2 

Finally, for resonant amplitude problems we will require 
C13 from 72-1 = C13723 + cuTZa. We see that in the near-field 
region M <C r <C cj -1 , 72. 1 is dominated by the growing- 
outward (r +2 ) solution while 72.3 and 724 are both domi- 
nated by the growing- inward (r 1- *) solution. Therefore the 
ratio C13 : C14 can be obtained by forcing the leading terms 
inward (i.e. coefficients of oc r 1 ' 1 ) to cancel. This is 

ki . fc 4 H .- e+1 (2£+ 1)! !_ t 

Cl3 = -^ Cl4 = 1 27^ = 1 } kl - (C 6) 



72a (r) 
72 B M 
72 c (r) 
72 D (r) 



iFi(^ + 3;2^ + 2;2io;r), 



e+2 - 
r e 

r 1 -V i " r ^{2- £;-2£;2iLor), 
r* +2 e lur 1 F 1 (£-l;2£ + 2;-2iLur), and 
r i-V" r (-2 - £■ -2£; -2iwr). 



(CIO) 



17 The product is empty for < 
to evaluate to unity. 



2, in which case it is understood 
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